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ABSTRACT 

We present models of ohmic heating in the interiors of hot jupiters in which we decouple the interior and 
the wind zone by replacing the wind zone with a boundary temperature 7; so and magnetic field B^q. Ohmic 
heating influences the contraction of gas giants in two ways: by direct heating within the convection zone, and 
by heating outside the convection zone which increases the effective insulation of the interior. We calculate 
these effects, and show that internal ohmic heating is only able to slow the contraction rate of a cooling gas 
giant once the planet reaches a critical value of internal entropy. We determine the age of the gas giant when 
ohmic heating becomes important as a function of mass, 7j so and induced B^q. With this survey of parameter 
space complete, we then adopt the wind zone scalings of Menou (2012) and calculate the expected evolution of 
gas giants with different levels of irradiation. We find that,with this prescription of magnetic drag, it is difficult 
to inflate massive planets or those with strong irradiation using ohmic heating, meaning that we are unable to 
account for many of the observed hot jupiter radii. This is in contrast to previous evolutionary models that 
assumed that a constant fraction of the irradiation is transformed into ohmic power. 

Subject headings: planets and satellites: magnetic fields — magnetohydrodynamics (MHD) — planets and 
satellites: atmospheres 



1. INTRODUCTION 

The large radii of many hot jupiters has been a puzzle ever 
since the discovery of the first transitin g planet HP 20 9458b 
( [Charbonneau et al.| (pOOO"); see |Baraffe et al.| ( [20T0l > for a 
review). |Guillot & Showman (2002) pointed out that if a cer- 
tain amount of the irradiation from the star (~ 1% of the inci- 
dent stellar flux) can be deposited deep in the envelope of the 
planet (pressures of > 10 bars), then the inflated radius can be 
explained. But the physical mechanism by which the required 
energy is transported into the interior of planet is still an open 
question. Several explanations have been proposed, such as 
a downward kinetic flux due to atmosphere circulation, or 



turbulent transpo rt and shock heating i n the flow (|Guillot & 
|Showman|[2002l [Youdin et al.|[20T0l JPern a et 
tidal dissipation (Bodenheimer et aT][2001 



of these has problems in accounting for al 



w (Guillot 



or 



2003 ). But each 
of the observed 



, 2011). 



r adii (e.g.|Laughlin e t al. 201 1 , |Demory et al. 

Batygin & Stevenson| (|2010) suggested that ohmic heating 
could serve as the heat source in the interior of inflated hot 
jupiters. The idea is that the shearing of the planetary mag- 
netic field by the wind driven in the outer layers of the planet 
by irradiation generates an induced current that flows inwards, 
dissipating energy by ohmic dissipation in deeper layers (see 
also Liu et al.|2008) . |Perna et aL] ( |2010a|b[ ) also pointed out 
the possible importance of magnetic drag on the dynamics of 
the flow in the envelope, and found that a significant amount 
of energy could be dissipated by ohmic heating at depths that 
c ould influence the rad ius evolution of the planet. 

|Batygin et al.| ( |201 1) implemented ohmic heating in evo- 
lutionary models of gas giants and showed that the amount 
of inflation depends significantly on the amount of irradiation 
received by the planet (and therefore its equilibrium tempera- 
ture r eq ). They found that the radii of low mass planets could 
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run away, increasing dramatically in response to ohmic heat- 
ing and leading to evaporation of the planet. This same be- 



havior was not observed in the recent study of Wu & Lithwick 
(2012), who find that ohmic heating could increase the radius 



of a planet that had already cooled, but only modestly. On the 
other hand, i f ohmic heating operates early in the lifetime of 
a hot jupiter, | Wu & L ithwick (2012) find that contraction can 
be halted and large radii obtained, large enough to explain the 
observed radii of all except a few planet s (see their Fig. 4 ). 

Both the time evolu tion models of Batygin et al. (20fTb and 
Wu & Lithwick ( 2012[ l calculate the profile of ohmic heating 



as J 1 1 a per unit volume inside the planet (where J is the cur- 
rent density and a the electrical conductivity), but adjust the 
overall level of heating so that the efficiency e — the frac- 
tion of the irradiation that goes into ohmic heating — is fixed 
at a level of e ~ 1%. In reality, the magnitude of the cur- 
rent that penetrates into the interior depends on how the flow 
in the wind zone interacts with the planet's magnetic field 
and the feedback from the magnetic field on the flow dynam- 
ics. |Menou| ( [2012| l considers scaling arguments for the atmo- 
spheric flows in a magnetized atmosphere, and argues that the 
efficiency e must decline at large T eq as magnet ic drag lim 



its the flow velocity in th e atmosphere (see also Perna et al. 
|2010bt [Rauscher et al.|2012) 



In this paper, we take a more general approach to the ques- 
tion of inflation due to ohmic heating, with t he aim of un- 
derstan din g the different ev olutions found by Batygir Tet al.| 
(2011 1 and Wu & Lithwick (2012), and incorporating the ef- 
fect of magnetic drag, and therefore variable efficiency, on the 
evolution. We take a different approach by separately con- 
sidering the planet interior and the wind zone. We first cal- 
culate the interior heating by replacing the wind zone with 
a boundary condition which specifies the toroidal field, or 
equivalently the radial current, at the base of the wind zone 
and the temperature there. This allows us to survey the pa- 
rameters that influence the amount of ohmic heating and its 
effect on the evolution of the planet. In this way we go beyond 
the previous assumptions of constant heating efficiency. We 
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then implement the scaling laws derived by Menou ( 2012| l, 
and show that indeed the efficacy of ohmic heating is reduced 
at high r eq because of increased drag in the wind zon e. In this 
way our tim e de pendent models differ c rucially from Batygin 
[eTaTlpOTTT l and |Wu & Limwick|p7T2] > in that we find that it 
is difficult to explain the observed radii of many hot jupiters 
with ohmic heating under the influence of magnetic drag, par- 
ticularly those with large masses M > Mj. 

The plan of the paper is as follows. In ^2j we review the 
general mechanism of ohmic heating and how the internal cur- 
rent is calculated, giving some order of magnitude estimates 
for the total ohmic power. Next in ^3] we present quasi- 
steady state models of gas giants undergoing ohmic heating 
as a function of their internal entropy S and the induced mag- 
netic field Z?0o and temperature 7i so at the base of the wind 
zone. We then use these models to follow the time evolu- 
tion of a given planet under the action of ohmic heating in ^4] 
With this general survey of pa rameter spa ce in hand, we then 
use the scalings derived by [Menou (2012! for the wind zone 
to calculate the evolution of observed hot jupiters and com- 
pare with observed radii (^5}. We discuss the limitations of 
our models and compare to other work in ^6] 

2. THE GENERAL MECHANISM OF OHMIC HEATING 

In this section, we review the basic physics of ohmic heat- 
ing, focussing on the generation of the induced field in the 
wind zon e an d radial current that penetrates into the planet 
interior (5 2.1i. We then estimate the likely magnetic field 
strengths that can be generated in the wind zone and the re- 
sulting ohmic power available for inflation ( §2. 2} . 

2.1. Calculation of the induced field and current distribution 

First we review the g eneral pict ure put forward by |Baty-| 
|gin & Stevenson] ( [2010) (See also ( jPerna et al |2010ab ). Due 
to strong irradiation from the host star, the hot jupiter has 
a thermally-ionized atmosphere, coupling the magnetic field 
and the atmospheric flow. The magnetic field could be either 
the intrinsic planetary magnetic field generated by a dynamo 
in the deep interior, or the external field from the host star. 
In either case, the evolution of the magnetic field in the wind 
zone is governed by the induction equation 



dB 



= -V x r/(V x B) + V x (v x B) 



(1) 



where v is the wind velocity and T] is the magnetic diffusiv- 
ity of the atmosphere. Assuming a steady state, and that the 
planet magnetic field in the outer layers is well-represented 

by a curl-free dipole field fidip, the induced magnetic field b is 
given by 

V x 7?(V x b) = V x (v x fi di p). (2) 

If the wind is predominately a zonal flow v* = v^cj), the in- 
duced field is toroidal, and will penetrate into the interior of 
the planet, with associated poloidal currents that close in the 
interior given by J = (c/47r)V x b. The internal toroidal field 
is given by 

V x (77V x = (3) 

below the wind zone where velocities are negligible. For a 
given poloidal current J, the local ohmic dissipation rate is 
^ohm = J 2 /cr, where a is the electrical conductivity, related to 
the magnetic diffusivity by 77 = c 2 /4na. 
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FIG. 1 . — The general mechanism of ohmic heating illustrated in a plane- 
parallel model. A vertical field B- is sheared by the wind v v in the wind zone. 
An induced field b is produced by the shear that penetrates into the interior. 
We use the temperature 7j so and induced magnetic field B^q at the base of the 
wind zone as boundary conditions for our interior solutions. 

Some intuition about the solution can be obtained by con- 
sidering a plane parallel model, which we illustrate in Fig- 
ure [T] There we divide the planet into three layers, repre- 
senting the outermost isothermal layer, with pressures lower 
than 60 mbars, the wind zone, between 60 mbars and 10 bars 
and the interior of the planet. The vertical field B z is sheared 
by the wind with velocity v v (z)exp(/fct). Focusing on the in- 
nermost layer representing the planet interior, and assuming 
constant conductivity, equation (Bl) gives an induced field b oc 
exp(-fe + ikx)y and associated vertical current 4ir j z /c = ikb. 
Both the field and current decrease exponentially in the inte- 
rior on a length scale 2ir / k. When the variation of conductiv- 
ity with depth is included, we must solve 



d 2 b 
dz 2 ' 



-k 2 b + d ^ 
dz dz 



= 0, 



(4) 



in which case the thickness of the penetration depth depends 
on the length scale on which the conductivity varies. 

This simple model illustrates that the interior solution de- 
pends on three factors: the geometry of the shearing veloc- 
ity in the wind zone (which is described by k in this simple 
model), the magnitude of the induced field at the base of the 
wind zone, and the profile of the electrical conductivity in the 
interior. The approach that we pursue in this paper is to con- 
sider the first two factors as boundary conditions on the inte- 
rior. We choose to parametrize our models by specifying the 
toroidal magnetic field at a pressure of 10 bars, which we will 
refer to as B^- 

To calculate the distribution of currents inside the planet in 
detail, we consider a simple geometry with a dipole field and 
a zonal flow v = vosinf?0. To solve equation (Bb in spherical 
coordinates we write the induced toroidal field: 



Ba, = sin # cos ( 



(5) 



in which case the radial dependence part of the field is given 
by 



V)- 



d\na , 



g(r) 



g'(r)-l(l + l) c ^=0 



(6) 



dr r 

where I = 2 is the index of associated Legendre polynomial 
. Having found g(r) and therefore B$(r), the currents are 
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determined by Ampere's law J = (c/47r)V x B. The ohmic 
power in the interior is 



'J 2 

P = I —dV 

a 



• 4tt 



{J? 



r dr, 



(7) 



where (J) = (J 2 +Jl) l l 2 is the effective angle averaged cur- 
rent at radius r. To get some feeling for the dependence of 
the field and current on depth, we can consider a power-law 
dependence a oc r a . In that case, the solution is oc r@ with 



(a-l)+\/(l+a) 2 + 24 



(8) 



and (J) 



,/8-l 



For more complex wind geometries, which 



involve I > 2, the solution is B^ oc r l Pj(cos6) for constant 
conductivity. For example, this indicates that more zonal jets 
in the wind zone implies a shallower penetration depth for the 
induced field. 

In fact, as we argue in the next section ( ^2.2[ ), the internal 
heating is dominated by the lowest densities, since the local 
heating rate is inversely proportional to the electrical conduc- 
tivity which increases rapidly with increasing pressure. This 
means that the current J can be taken as a constant without 
making a significant error in the heating profile. We have 
confirmed this by comparing constant current solutions with 
detailed solutions of equation (|6j. 

For the constant current case, we compute the current as 



4^ ^7 



(9) 



where Rj = 7 x 10 9 cm is the radius of Jupiter, B^q = B^i^p = 
lObars). This means that there is a direct mapping between 
the chosen value of B^o, the radial current inside the planet 
J r , and the local heating rate, taken to be J 2 /a(r) per unit 
volume. Note that we do not take into account the averages 
over angle in equation dTli nor the true radius of the planet, 
and so our value of B^q for a given amount of internal heating 
could be a factor of a few of the toroidal field in a model 
which self-consistently includes both the wind zone and the 
interior. Instead, our parameter B<j& should be interpreted as 
a measure of the internal heating (given by eqs. |9) and then 
J 2 . 1 a locally). 

2.2. Magnitude of the induced field and ohmic power 

It is useful to estimate the expected magnitude of ohmic 
heating and how it scales with parameters such as planet mass 
M. The first step is to use equation Q to estimate the expected 
strength of the induced field by dimensional analysis, 



b = B, 



dip 



a 

lOV 



4naHv 



(10) 



OMRj) (lkms- 1 ) ' (U) 



where v is an average wind speed and a a typical value of 
electrical conductivity in the layer] and we take the vertical 
length scale to be the pressure scale height H. 

3 We give the electrical conductivity in cgs units here. Note that the con- 
version to SI units is 1 S irr' = 9 X 10 9 s _I . 



In this paper, we take Z?dip = 10G as a standard value. Typ- 
ical dipole field strengths for hot jupiters have been esti- 
mated from scalings with planet parameters (see Trammell 
|et al.|20fT|fo r a detailed summary and discussion). Sanchez- 
|Lavega| ( |2004| ) argued that the field is generated by the dy- 
namo action in the metallic region as in Jupiter ( jStevenson] 
|1983| l, with the field strength closely related to the rotation 

of the planet, B ~ (pflrj)^. This predicts that the field on 
typical hot jupiters should be a factor of few smaller than 
that on the Jupiter, with a typical value of equa torial mag- 
netic field B eq ~ 5G. However, Christ ensen et al.| ( |2009[ ) ar- 
gue that the field instead scales with the heat flux escaping 

from the conductive core at large enough rotation rate, giving 

9 i 

B ~ (p,F C o re )5. This gives a field strength an order of magni- 
tude larger than estimated with the previous method. 

Given an induced field b, we can then estimate the expected 
ohmic power per unit mass, P m = (J) 2 /(pa). Approximating 
the angle averaged current at the base of the wind zone using 
the constant current case as equation Q gives 



P m =10- 1 ergg" 1 s 



B<po V 

Tog/ 



io 6 s- 



10" 4 g cm" 



(12) 

Since a increases dramatically in the core of the planet, heat- 
ing at low density dominates. If the conductivity profile scales 
as exp(-r/H) for example, where we take the lengthscale as 
the pressure scale height H, the total power in the interior is 



^ohm. total 

= 10 23 erg s" 1 



(13) 
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(14) 



where the subscript t indicates a quantity evaluated in the out- 
ermost regions of the interior just below the wind zone. Note 
that H = KT/g, where we adopt K = 3.64 x 10 7 erg g _I K _I 
for a hydrogen molecule dominated composition with helium 
fraction Y = 0.25. As we mentioned previously, because the 
total ohmic power is dominated by the heating at low pres- 
sure, it is not very sensitive to the radial profile of the current. 
The scaling in equation ( 13 1 implies that P h m oc 1 jM when 



the radius and conductivity of planet is not strongly depend 
on mass, a scaling that we find in our numerical solutions. 
More massive planets have less ohmic power for a given B^q 
and temperature 7; so at the base of the wind zone. 

3. OHMIC HEATING AS A FUNCTION OF INTERNAL 
ENTROPY 

In this section, we calculate the structure, luminosity, and 
ohmic heating profile for gas giants as a function of their in- 
ternal entropy S. We will use these models in 34] to follow 
the time evolution of the planet by following the aecreasing 
entropy as the p lanet co ols. In this approach, des cri bed by 
Hubbard ( 1977 ) (see also |Fortney & Hubbard|2003| and [Arras 
& Bildsten||2006[ ), it is assumed that the convective turnover 
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time is much shorter than the evolution time of the planet, so 
that the convection zone maintains an adiabatic profile as it 
cools and lowers its entropy. We also assume that the radia- 
tive envelope has a thermal timescale much shorter than the 
evolution time, so that the envelope is in thermal steady-state, 
carrying the luminosity emerging from the convection zone. 
The luminosity of the planet is then given by the radiative lu- 
minosity at the radiative-convective boundary, 



L = 



l6irGcM r aT 4 
3np c 



(15) 



where p c is the pressure at the convective boundary, T is the 
temperature at that location, and M r is the enclosed mass. The 
evolution of the internal entropy is then given by 



dS , 
I — dm ■ 
dt 



MT- 



,dS 
dt 



-US) 



(16) 



where we have assumed dS/dt is constant across the con- 
vection zone and define the mass-averaged temperature f = 
J T/Mdm. 

The effect of ohmic heating appears in two places in this 
approach. The first is that an ohmic heating term must be 
added to the right hand side of equation ( [To} , 

J 2 



MT- 



,dS 
dt 



= -L r + 



dV, 



(17) 



where the integral is over the convection zone. For a planet 
with a given entropy S, the luminosity at the top of the con- 
vection zone is fixed by the structure and is given by equation 
(fl5). However, because some of this luminosity is now pro- 
vided by ohmic heating, the cooling rate of the convection 
zone (dS/dt) is smaller. The second influence of ohmic heat- 
ing is that it can change the temperature profile in the radiative 
zone, in particular by pushing the radiative-convective bound- 
ary to higher pressure and lowering the luminosity (eq. fT3j|). 

In this section, we include the first effect by calculating 
gas giant models without feedback from ohmic heating in 
the atmosphere (§3.1), and then include the feedback from 
ohmic heating in the radiative zone to include the second ef- 
fect (§3.2). In §3.3, we summarize the results. 

3.1. Planet models without feedback 

We now make models of gas giants with given mass M and 
central entropy S, and use them to calculate the ohmic dissi- 
pation in the planet, but without including the effect of ohmic 
heating on the planet structure. This allows us to calculate the 
ohmic heating within the convection zone and, by including 
this ohmic power in equation ( 16 1, the corresponding slowing 
of the cooling rate. 

We present the detail microphysics in our planet model 
in the Appendix [A] Here we only note the differences be- 
tween our opacity and conductivity profiles and those used 
in other works. Our opacity profile use a different extrap- 
olation method in the intermediate pre ssure range between 



10 5 bars from Paxton et al. 



(2011 1, who take the opac- 



10 3 - 

ity for log R > 8 (where R = p/T£ is used in the opacity 
tables) to be a constant set by the value at log/? = 8. As 
far as we are aware, opacity calculations for this interme- 
diate pressure range have not been carried out. For plan- 
ets undergoing a large amount of internal heating, the con- 
vection zone boundary moves into this region, and so know- 
ing the opacity there is important for understanding the loca- 
tion of the convection zone boundary. As we describe later, 
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FIG. 2. — The temperature profile (top panel) and entropy profile (bot- 
tom panel) for different treatments of the outer layers, either isothermal (red 
curve), radiative (green curve), or radiative including ohmic heating (blue 
curve). For the ohmic heating model, we take B^q = 1000 G. In each case, 
the radiative-convective boundary is marked with a vertical bar. Note that we 
do not show the entire structure, but focus on the lower pressures to illustrate 
the differences in the position of the radiative-convective boundary between 
models. 



this controls the contraction rate of the cooli ng planet. Our 
potassiu m conductivity prof ile (equation |A2| ) is t he same as 
used by|Perna et aT] ( |2010a| , but is different from |Batygin &| 

n 

rather than 10~ 15 cm 2 



iStevensonj ( |2010[ ), who use a electron-neutron cross-section 
of tt(7.2 x lO^cm) 2 = 1 .6 x 1(T 16 1 
( Draine et al.||1983) and a slightly different thermal averag- 
ing factor for th e velocity. The resulting diffe rence is that the 
conductivity of Batygin & Stevenson (2010) is a factor of 9 
times larger than our conductivity. 

The planet model is calculated by integrating outwards 
from the center, following the convective adiabat until it in- 
tersects a radiative zone extending inwards from the surface. 
When calculating the cooling curve of an irradiated gas gi- 
ant, a reasonable approximation is to take the outer radiative 
zone of the planet to be isothermal. However, because ohmic 
heating is very sensitive to pressure, a small error in the deter- 
mination of the convective boundary results in a much larger 
error in the ohmic power. To illustrate this, we have calculated 
models with an isothermal radiative layer and with a radiative 
layer that is in thermal equilibrium and carries a constant lu- 
minosity equal to the luminosity from the convection zone. 

For the isothermal case, we integrate 



dm 



= 4npr 



(18) 
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TABLE 1 
Model Summary 



Model 


R(Rj) 




Pconv (bars) 


^ohra(erg S ')(/)> pconv) 


Pohm(erg s ')(p > 10 bars) 


iconv(erg s ') 


Isothermal 


1.25 


10 


62.8 


8.0 X 10 23 


1.2 x 10 25 


1.3 x 10 26 


Isothermal 


1.25 


100 


62.8 


8.0 X 10 25 


1.2 x 10 27 


1.3 x 10 26 


Radiative 


1.25 


10 


131.7 


6.8 x 10 22 


2.8 X 10 24 


7.7 x 10 25 


Radiative 


1.25 


100 


131.7 


6.8 x 10 24 


2.8 x 10 26 


7.7 x 10 25 


RadiativeFB b 


1.25 


10 


132.0 


6.7 x 10 22 


2.8 x 10 24 


7.7 x 10 25 


RadiativeFB b 


1.25 


100 


176.4 


3.1 x 10 24 


1.9 x 10 26 


5.6 x 10 25 



a Model computed with 5 = 8, 7i s0 = 1500 K, and M = 0.96 Mj. We refer to this set of input parameters our standard model 
in the text. 

b Model including the feedback of ohmic heating in the atmosphere. 

c Both the cooling luminosity and the ohmic heating in the convective zone reduce while the convective zone boundary 
move towards deeper pressure due to the feedback in up atmosphere. 
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FIG . 3 . — The cumul ative ohmic power against pressure for the same planet 
parameters as in Figure fTT] with B = 10G (radiative model, no feedback). The 
solid curve uses a current profile calculated by solving equation |6|; the red 
dashed curve assumes a constant current with depth. 



dp 
dr 

dT _T 
dr 



Gm 

r l 



P dr' 



(19) 
(20) 



outwards from the center, taking V = V a d for T > 7] so (adi- 
abatic interior) and V = for T < 7j so (isothermal layer). In 
the non-isothermal case, we take 



V rad : 



; min(V rad ,V ad ) 
3kL p 
\6-KcGM~aT*' 



(21) 
(22) 



We calculate a model for a given M and S by integrating 
outwards from the center and inwards from the surface to 
a matching pressure p = 30 kbars. For the outwards inte- 
gration, we choose the central pressure p c and cooling rate 
dS/dt (or equivalently cooling time t$ = S/ \dS/dt\). For the 
inwards integration, we start at a pressure of 10 bars and set 
the temperature there to be 7; so . We then integrate inwards, 
choosing the luminosity L and radius R. A multi-dimensional 
Newton-Raphson method is used to find the correct choices 
of (p c , ts,R,L) that result in m, r, T and L agreeing to within 
1% at the matching pressure. 
As an example, Figure [2] compares the entropy and tem- 



perature profiles for models with an isothermal and non- 
isothermal radiative zone. In the isothermal case, we choose 
T c = 3 x 10 4 K, p c = 2 x 10 7 bars and r iso = 1500 K, which 
gives a M = 0.96 Mj, R = 1.25 Rj planet with core entropy 
S = 7.98 and convective zone boundary p conv = 62.76 bars. 
The luminosity from the interior is L = 1.28 x 10 26 erg s" 1 , 
giving t s = 5.78 Gyr. With a radiative zone, we obtain the 
same mass, radius and entropy with a convective zone bound- 
ary Pconv = 131.7 bars. The luminosity from the interior is 
L = 7.7 x 10 25 erg s" 1 , and t s = 10.5 Gyr to cool. 

In each case, the magnetic field structure in the planet in- 
terior is obtained by solving equation (|6]l using the conduc- 
tivity profile in the planet (shown in Figure [17), and then the 
ohmic heating profile is determined. For an induced magnetic 
field Z?0 ( ) = 10G at the bottom of the wind zone (p = 10 bars), 
we find P hm(p ^ Pconv) = 8 x 10 23 erg s" 1 for the isothermal 
model and P h m (p ^ Pconv) = 6.8 x 10 22 erg s" 1 in the non- 
isothermal case. While the cooling time for the planet changes 
by a factor of two between the two models, the ohmic power 
changes by more than an order of magnitude. Therefore, it 
is crucial to locate the convective boundary accurately when 
calculating the ohmic power inside the convection zone. 

In Figure[3] we compare the ohmic power calculated in this 
way, which includes the correct radial distribution of current, 
with the ohmic power calculated by assuming a constant ra- 
dial current, independent of depth. The agreement is excel- 
lent (within a factor of 2) except at the highest pressures in 
the central regions of the planet. 

3.2. Planet models with feedback from ohmic heating 

The fact that the ohmic heating per unit mass rises rapidly 
to lower densities (Fig. [3]) suggests that the heating in the re- 
gions lying between the wind zone and the convection zone 
boundary will be larger than the heating in the convective in- 
terior. We include the ohmic heating in the radiative layer by 
allowing L to vary throughout the radiative zone, with 



dh _ J 1 
dr a 



(23) 



We do not include ohmic heating at pressures less than 10 
bars. Instead, we specify the temperature 7; so at p = 10 bars 
and integrate inwards. Of course, there could be significant 
ohmic heating within the wind zone at p < 10 bars, but we 
absorb this into the boundary condition. Note that this means 
that early in the lifetime of the planet, when the entropy is 
large enough that p conv < 10 bars, the models here revert back 
to our previous models with no feedback. However, at those 
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FIG. 4. — The luminosity L (solid curves) and the total ohmic heating in 
the convection zone P hm (dashed curves) as a function of central entropy S. 
The red, yellow, green and blue lines (from top to bottom for the solid lines; 
inverse for the dashed lines) represent planets with different mass: 3.0 Mj, 
1.0 Mj, 0.6 Mj, 0.3 Mj. At larger entropy, where ohmic heating is unimpor- 
tant, L oc M at fixed S, whereas P h ra decreases with increasing M. All the 
planet models are computed with 7j so = 1750 K and B^q = 100 G. 
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FIG. 6. — Position of the radiative-convective boundary as a function of 
internal entropy S for a 1 Mj planet with 7j so = 1750 K. Blue, green, yellow, 
red, black lines are for B^ =0, 10, 30, 100, 300, 1000 G. At a given entropy, 
a larger B^q results in more ohmic heating in the radiative zone, moving the 
convective boundary to a higher pressure. 




FIG. 5. — The time history of planet luminosity (solid curves) and ohmic 
heating in the convection zone P h m (dashed curves) for M = 0.3 Mj (blue), 
0.6 Mj (green), 1 Mj (yellow) and 3 Mj (red curves) (same configuration with 
Figurep). The luminosity decreases with time because of cooling, while the 
ohmic Seating either increases or decreases slowly depending on M. At late 
times, when ohmic heating in the radiative layer becomes important, P h m 
decreases because the convective boundary moves inwards. When P h m be- 
comes comparable to L, the cooling and contraction of the planet is halted. 
All the models are calculated with 7j so = 1750 K and B^, o = 10" G. 

early times, ohmic heating is generally not yet important. 
Also note that in the models without feedback, the strength 
of the induced field B^ does not influence the internal struc- 
ture of the planet, whereas here a larger B^ results in more 
heating in the radiative layer which can push the convective 
boundary deeper. 

In Table [1] we compare the models with feedback to our 
earlier models without feedback. The internal structures are 
shown in Figure [2] The planet radius does not vary much be- 
tween different models. The biggest difference is in the posi- 
tion of convective zone boundaries (marked by black vertical 



bars in Fig. Et, which results in a difference in the cooling lu- 
minosity ana the ohmic heating both in the convective zone 
and atmosphere. Note that this means that the cooling history 
of a planet using these three approaches would be different, 
especially the time at which ohmic heating begins to become 
important for evolution. We use this feedback model for all 
the calculation carried on below. 

3.3. Ohmic power as a function of entropy 

The luminosity and ohmic power is shown as a function 
of entropy in Figure |4] As noted in particular by Arras & 
Bildsten (2006), equation ( p"5j ) shows that L oc M at fixed en- 
tropy, and we see that scaling in Figure [4] On the other hand, 
the ohmic power decreases with increasing M, as discussed 
in §2.2. We find that the decrease in P h m is well described 
by fohm oc R 2A /M, which has a shallower dependence on R 
than in equation (|T3J because of the dependence of the con- 
ductivity term on mass which compensates the R 4 term. In 
our feedback model, for a lower mass planet the higher atmo- 
spheric ohmic heating pushes the convective zone boundary 
slightly deeper, resulting in a higher conductivity at the top 
of the convection zone. Overall, lower mass planets generally 
have higher ohmic power deposited in the convection zone. 

Combined with the mass-luminosity dependence, the de- 
crease of ohmic power with mass means that ohmic heating 
becomes important for lower mass planets at a much higher 
entropy than for more massive planets. The value of entropy 
at which ohmic heating becomes important depends on the 
boundary induced field B^q. Turning this around, for an ob- 
served planet with measured radius and mass, we can infer 
the entropy and therefore derive a limit on the wind zone B^o 
required for ohmic heating to be providing a significant part 
of the luminosity in that object. We carry out this procedure 
in ^5] but first describe our calculations of the time-evolution 
of planets with ohmic heating. 

4. TIME DEPENDENT EVOLUTION OF PLANET 
STRUCTURE 
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FIG . 7 . — The age of the planet when ohmic heating becomes important 
(iota > 0. IL) versus planet mass, for B 0O = 30 G (black), 100 G (red), 300 G 
(green), 1000 G (blue curves), and for two temperatures 7i so = 1750 K (solid 
curves) and 2250 K (dashed curves, corresponding to the lower panel of la- 
bels). 
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FIG. 8. — The value of B^q required for ohmic heating to become important 
at an age of 3 Gyr. From bottom to top, M =0.3, 1, 3.0 Mj. 



Having computed the luminosity at the radiative- 
convective boundary for a large grid of models with differ- 
ent M, S, T[ so and B^o, the evolution in time of a planet with 
fixed mass M, Bjq and 7| s0 



then involves stepping in entropy 
using equation (1171). When calculating the ohmic power, we 
assume constant current J with depth, which as shown earlier 
is a good approximation. We have checked that for B^q = 0, 
our cooling models c omp are well with the earlier results of 
Burrows et aL] ( [19971 ) and |Baraffe eTaT] ( [2003l > (see IMarleau] 
et al.|2012 1. We reproduce their cooling curves to within 30% 
in luminosity, and predict radii that are 0.05-0.08 Rj larger 
than in those cooling sequences. 

By integrating in time, we compute the time history of the 
planet luminosity and ohmic power, shown in Figure [5] As 
the planet cools, the convection zone ohmic power increases 
or decreases slightly depending on mass, but always changes 
more slowly than L, so that ohmic power eventually becomes 
comparable to the cooling luminosity. At the same time, as 



FIG. 9. — Time history of planet radius for (top to bottom) M =0.3, 0.6, 1 
and 3 Mj. All the models are calculated with 7; so = 1750 K and B^o = 100 G. 



the ohmic heating in the upper atmosphere (which is about an 
order of magnitude larger than convection zone heating) starts 
to affect the planet structure, the convective zone boundary 
shrinks inwards. The resulting decrease of both cooling lu- 
minosity and convection zone ohmic heating result in a rapid 
increase in cooling time, so that we can view the evolution 
afterwards as a quasi-steady state. For higher atmospheric 
ohmic heating, this effect happens at higher entropy, thus the 
steady radius of planet will be larger. We compare the evo- 
lutionary tracks of the radiative/convection zone boundary of 
a 1 Mj planet with different strengths of induced field in Fig- 
ure [6] As we increase the amount of ohmic heating, the con- 
vection zone boundary deviates from the no heating path at a 
higher entropy. 

In Figure [7] we show the age of a planet with a particular 
mass when P h m ,t' = 0.1L conv , at which point ohmic heating 
starts to become significant and the planet contraction slows. 
In general, the atmospheric heating is an order of magnitude 
larger, and so comparable to the cooling luminosity at this age. 
In FigurelTJ we report that the result is sensitive to the strength 
of B^. With r iso = 1750 K, for B^ equals 30 G, a 0.3 My 
hot jupiter can reach steady state in 1 Gyr, while a 1 Mj 
requires 100 G to reach steady state at a similar age. Higher 
7] so (dashed line in Figure |7| does not help the planet reach 
a steady-state radius faster, out on the contrary, it requires a 
larger induced field to achieve the same result. 

Figure [8] shows the B^o required to halt contraction within 
3 Gyr as a function of 7; so . We see that a hotter planet re- 
quires a stronger induced field to obtain a significant level of 
ohmic heating. This is because the interior ohmic heating is 
closely related to the conductivity at the bottom of the wind 
zone; the conductivity increases with temperature, reducing 
the ohmic heating at fixed induced field. But we should also 
point out that for a hotter planet there is a higher chance to 
obtain a stronger induced field due to the stronger wind in the 
atmosphere. So this result does not necessarily imply that it 
is more difficult to make ohmic heating important in hotter 
planets. 

We plot the time evolution of the planet radius for plan- 
ets with B^o = 100 G and 7j so = 1750K and different planet 
masses in Figure [9] In the absence of stellar ages, we shall 
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take a typical age of 3 Gyr, and take the radius at 3 Gyr as 
the present day radius. For the lowest mass planet, 0.3My, 
the effect on the evolution of the planet radius is significant, 
and the planet stops cooling around 1 Gyr (with cooling time 
longer than 10 Gyr) and thereafter maintains a large radius. 
However, the heating is not as effective at higher masses. For 
example, HD 209458b has a observed radius of 1.35 Rj with 
mass 0.7 My, while we can only obtain a radius of 1.25 Rj 
for ZJ^o = 100 G. This is because the power we introduced 
into the planet interior is far smaller than the received stellar 
luminosity. In the case of our standard model, the irradiation 
luminosity from the host star is 10 29 erg s" 1 , and the heating in 
the interior is only one 0.01% of it, 10 26 erg s -1 . To go further, 
we must understand what values of B^q might be expected as 
a function of T eq , and we turn to this in the next section. 

5. EVOLUTION INCLUDING WIND ZONE MODEL 
AND COMPARISON TO OBSERVATIONS 

In ^4] we calculated the time-evolution of cooling gas gi- 
ants assuming that 7i so and B^q are independent parameters. 
In reality, they are coupled by the dynamics in the wind zone, 
since the atmospheric flow, in response to the irradiation, de- 
termines both the magnetic field in the layer and the tem- 
perature at depth (the values of 7j so and B^o are specified at 
p = 10 bars). In this section, we imple ment th e scaling s for 
the wind zone dynamics proposed by Menou (201 2) (§5.1) 
and then compare our results to observed systems ( §5. 2} . 

5.1. Dynamics of the wind zone and the relation between 
7j so and Bm 



Both |Batygin et al.| ( |201 1) and |Menou"| ( |2012J l write down 
simplified models for the wind zone dynamics i ncluding the 
effects o f magnetic drag. In both cases, following Perna et al. 
( 20 10a} the magnetic drag force is assumed to be J x B/c per 
unit volume, with the current J set by a balance between the 
shearing of the magnetic field by the fluid and ohmic diffusion 
of magnetic field lines against the fluid motion, J = crv x B/c 
(^2|. However, the dynamical balance in the two models is 
quite different. Menou (2012 1 writes the force bal ance for the 
equatorial flow as (see also Show man et al.|2 010) 



= - 



Rf 



UAT hmiz A\np v^B 2 
Rp Aitprj 



(24) 



The first two terms represent a balance between the advective 
term and the horizontal driving from the day-night tempera- 
ture difference A7h or i z . This balance is thermal wind-like in 
that the horizontal pressure gradients require a vertical gradi- 
ent in the fluid velocity over a vertical pressure scale A In p. 
The final term represents the magne tic drag force, again i nte- 
grated over a vertical scale A In/?. Batygin ~et al.| ( |20TT) on 
the other hand consider the meridional circulation induced by 
magnetic drag on the azimuthal flow, so that for example the 
latitudinal force balance is fv y = /tt_ where / = 2ft sin 8 is 
the Coriolis parameter and n the magnetic drag timescale. 
Their solution represents a thermal wind balance involving 
the equator-pole temperature gradient, modified by magnetic 
drag. 

In both cases, magnetic drag limits the fluid velocity at 
high temperatures, where the large degree of ionization and 
therefore large electrical conductivity results in strong cou- 
pling of the fluid and magnetic field. Balancing the second 
and third terms in equation ( 24 1 gives oc 77 when magnetic 



drag dominates, and therefore the magnetic Reynolds number 
Rm = v^H/r) becomes almost constant, varying only slowly 
with temperature. Similarly, equation (16) of Batygin et al. 



(201 1[) has two possible limits, either oc 77 when the lat 



eral temperature gradient is large, in which case R M becomes 
almost constant at large 7i so , or oc i] 2 when the drag time 
scale is comparable to the rotation period, while the lateral 
gradient of temperature is still small, in which case Rm oc rj 
declines rapidly at large 7j so . 

A dynamical model including both day-night driving and 
meridional circulation with magnetic drag is not yet available. 
For ou r purposes, we have implemented the model of Menou 
(20121 as described by equation ( |24} , with the day-night tem- 
perature difference given by 



(5 s ;) T adv < r rad 



(25) 



In equation ( 25 1, 7d ay is the dayside averaged temperature 



considering a dilution factor of 0.5, = 2T^ y = 47^, and 
the advective and radiative timescales are 



7"adv ' 



7"rad : 



Rp 

C P p 
8^sbT 6 3 



(26) 



(27) 



day 



Note that these timesc ales are evalu ated at the outermost pres- 
sure, which following Menou (2012) is taken to be 60 mbars, 
the estimated location of the thermal photosphere. This is 
the reason for adopting the thermal timescale appropriate for 
an optically thin region, so that r ra( j oc p in equation ( |2"7] >; in 
deeper, optically thick layers, r, al j has an extra factor j)f the 
leading to r ra( i oc p 2 



optical depth 



ma n et al 
60 mbars 



(e.g. Fig. 3 of |Show- 
The magnetic drag term is also evaluated at 



p = ou moars; this term is integrated over height, but since a 
decreases with increasing pressure (for an isothermal layer), 
the dominant contribution to the integral is from the lower 
limit on pressure, and so 77 and p ar e evaluated th ere. This is 
an important difference from Baty gin et al.| ( |201 T) , who evalu- 
ated their magnetic drag timescale at p = 10 bars, which gives 
a drag time an order of magni tude longer tha n we find here. 
Based on that estimate, Batyg lrTet al.| ( 20TT) concluded that 
the drag timescale was always much longer than a rotation 
period. 



We solve equation (24 1 for as a function of 7" eq , and find 
the corresponding value of the induced field B^ from equation 
( 10 1 for different values of the dip ole field ffd ip- For A In/? = 
O.y, we reproduce the results of Menou (2012) (see his Fig. 1), 
but we also consider larger values of A In/?. |Menou| ( |2012| l 
models the weather layer with a modest vertical extension 
around 1 pressure scale height. We also solve the equation 
with A In/? = 3 for typical values in hot jupiter atmosphere as 
reported by |Showman et al" ( 2010 1, and Aln/? = 5 fo r a wind 
zone extending to p ~ 1 bars, for compariso n with Batygin 
|& Stevenson| ( p0T0l i and [Batygin "etaLl ( pOTT) . The effect of 
varying Aln/? on B^ is shown in the left panel of Figure 10 
For numerical convenience, we fit the B, 
following: 



T eq relation with the 



1 _ 1 1 

B^,(T eq ) B a( jv -Bdrag 



(28) 
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FIG. 10. — The induced field (left panel) and timescales (right panel) in the wind zone as a function of T e q. In the left panel, the solid, dashed and dotted 
curves are for wind zone thickness Alnp = 0.9,3 and 5, and we take B r = 10 G. In the right panel, we show the advection, radiative and drag timescales ijj,, 
r ra( ] and r^ao. The advective timescale is shown for A lnp = 0.9 (solid), 3 (dashed) and 5 (dotted curves) (the radiative and drag timescales are independent of 
Alnp). 



where 



5adv = 2.8x 10 6 Gr eq exp 



2.53 x 10 4 



m" 2 



B r 

10 G 



and 



B 



drag 



= 1125 G 



(-. 



V 1000 K 



Alnp 



10 G 



(29) 



(30) 



Changing the wind 
5 moves the tran- 



This reproduces to within ss 10% for r eq in the range 1 100 
to 2200 K. 

The transition to the regime where Rm is approximately 
constant occurs when Td rag exceeds r a d v , where the magnetic 
drag timescale is Td rag = 4-np?]/B 2 r = rj/v\, where va is the 
Alfven speed, and again the timescale is evaluated at the top 
of the wind zone (p = 60 mbars here). These timescales are 
plotted as a function of T eq in Figure 10 
zone thickness from Amp = 0.9 to Amp 
sition temperature from r eq s» 1400 K to 1700 K. 

To use this value of as a boundary condition for our evo- 
lutionary models, we must relate the temperature T eq at low 
pressure to the temperature 7| so at p = 10 bars. This relation 
depends on the details of energy transport in the wind zone, 
including the effects of ohmic heating and needs to be studied 
furthe r. Here, we a dopt the atmospheric temperature profile 
from Guillot (2010), and keep in mind the uncertainty in the 
relation between r eq a nd 7j so when i nterpreting our results be- 
low. The relation from |Guillot] ( |20T0| is (see his eq. [29]) 



37^ 



2 1 
3 + 



(31) 



where 7 is the ratio between visible and infra-red opacities, 
and f = 1/2 for a dayside average or / = 1/4 for an aver- 
age over the whole surface. Choosing 7 = 0.4, as appropriate 
for a p lanet like HD 409658b (e.g. see Fig. 1 of Hubeny et al. 
2003) and a dayside average /= 0.5, we obtain T[ so = 0.947j rr = 
1.33r eq . In the following section, we will use this relation to 
infer the appropriate value of 7[ so from the T eq of observed 
planets. We note here that we don't have a good knowledge 
of what the 7 parameter would be for most of the observed 
planets. While 7 parameter could vary in a very large numeri- 
cal range, the ratio between 7] so and T eq only changes within a 
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FIG. 1 1 . — The predicted mass-radius relation at 3 Gyr for B^q =0, 30, 100 
and 1000 G and T eq = 1316 K (solid curves) and 1692 K (dashed curves). The 
data points show observed transiting planets, divided into two temperature 
groups T > 1500 K (green points) and T < 1500 K (red points). 

factor of few [0.99(7 ~> inf) < CW r eq) < 3.05(7 = 0.01)]. 
Since the observed properties of planets gives T eq thus the 
boundary induced field, changing the value T\ so /T eq is equiva- 
lent to shift inside the plane T^-B^ given by Figure [8] Gen- 
erally, a smaller 7i so /r eq is favored to inflate the planet with 
the same B^. 

5.2. Comparison with Observed Hot Jupiters 

In Figures [TTJ to JT3] we compare our results with the ob- 
served properties oftransiting planets taken from the TEPcat 
transiting planet catalog which gives the planet mass, radius, 
and equilibrium temperature r eq = 7^, e ff {R*/2a) l l 2 where T^eff 
is the stellar effective temperature. As the ages of most stars 
are unknown or highly uncertain, we assume an age of 3 Gyr 
when comparing with the observed planets. 



First, Figure 1 1 shows the effect of increasing Bm at fixed 



7j SO on the planet radius. To help compare with the data, 
we divide the observed sample into two groups with either 
r eq > 1500 K (green points) or < 1500 K (red points) and 
show theoretical curves for either T eq = 1316 K or 1692 K 
(these two temperatures correspond to T lso = 1750 and 2250 K 
respectively). We see that for the low r eq group, an induced 
field of 10-100 G can explain most of the observed radii, 

4 http://www.astro.keele.ac.uk/~jkt/tepcat/ 
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FIG. 12. — Comparison with observations using B^oC^eq) from the wind 
zone model. Top panel: Radius at 3 Gyr against r eq for M = 0.3 (red), 0.6 
(green), 1.0 (blue) and 3.0 Mj (black). The data points are observed planets 
divided by mass: 0.2 Mj < M < 0.5 Mj (red points), 0.5 Mj <M < 0.9 Mj 
(green points), 0.9 Mj < M < 1.3 Mj (blue points). Bottom panel: Predicted 
mass-radius relation at 3 Gyr for T e q = 1316 K and T = 1682 K (bottom to 
top). In each case, the dashed curve shows the radius without ohmic heating; 
the solid curve with ohmic heating. The data has been divided by tempera- 
ture: T eq < 1500 K (red points), T > 1500 K (green points). 



while the high T eq planets need at least B^q = 1000 G to match 
the observed radii. It is clear that a higher induced magnetic 
field is needed to explain a given radius at higher equilibrium 
temperature. 

Next, we use the wind zone model described in §5.1| to cal- 
culate B^o as a function of T eq , assuming canonical values 



B, = 10 G and Alnp = 3. In the top panel of Figure 12 
show the radius as a function of T eq , with the colors represent- 
ing three different bins in planet mass. There exists a clear 
correlation between the radius and r eq , both in the observa- 
tions and the models. In addition, we see that the amount of 
inflation is also strongly dependent on the planet mass. Plan- 
ets within the mass bin 0.3-0.6 Mj agree quite well with our 
ohmic heating model. However, ohmic heating clearly can- 
not explain planets with mass ~ 1 My and large inflated radii 
> 1.4 Rj. Ohmic heating can help to increase the radius (for 
comparison the dashed line shows models with no ohmic heat- 
ing), but not enough to match the observed value. This is a 
consequence of the increased power needed to maintain the 
radius of a massive planet at a particular value, as well as the 
reduced ohmic heating power at larger masses. 
In the lower panel of Figure [T2j we show the radius as a 



function of mass. As in Figure [TT] we divide the data into 
two temperature ranges and showlne model results for two 
representative temperatures, now using the wind zone model 
to specify the value of B^q for each temperature. The param- 
eter region where ohmic heating has the largest effect is high 
temperature, low mass planets. The radii of the low tempera- 
ture group (T eq < 1500 K) can almost all be explained without 
ohmic heating. For the high temperature group, ohmic heat- 
ing can explain the observed radii of low mass planets, but 
most of the radii of the high temperature group lie well above 
the models, especially at large planet masses > 1 Mj. 

In Figure [13] we show the ratio between observed and pre- 
dicted radii R~a>s/Rprei against r eq (upper panel) and against 
M (lower panel) for each observed planet. In this case, we 
use the observed values of M and T eq to calculate the evolu- 
tion of the planet, and, in the absence of stellar ages, we take 
Rp le d to be the radius at 3 Gyr. In the upper panel, we see 
that for T eq > 1600 K, there are many planets whose radii lie 
above the predicted values. The slow increase of B^q with 71 s0 
at large 7; so due to the magnetic drag term results in a much 
weaker dependence of /? pre d on T eq than observed, and most 
outliers lie at the highest temperatures. In the lower panel, we 
see that the majority of the unexplained objects (R b s > Rpred) 
are at larger masses M > 0.7 Mj. We note that the choice 
of estimating the planet radius at 3 Gyr is not critical for the 
above picture. Constraining ourself within the time range of 
1-5 Gyr, the predicted radius only varies within several per- 
cent. 

To look in more detail at the effect of our assumed wind 
zone model on how successfully we are able to reproduce the 
observed radii, in Figure 14 we show the results in the B^q- 



T eq parameter space. For each observed planet, we first make 
a model with no ohmic heating, varying the internal entropy S 
at the measured M and T eq until we match the measured radius 
R p . Then we calculate the value of B^q required in that model 
for the ohmic power in the convection zone P ohm to be 30% 
of the planet's luminosity. We colour-code the data points 
according to whether they are successfully explained by our 
time evolutions, ie. whether they have /? pre d larger or smaller 
than /? b s in Figure[T3] These two groups of data points lie on 
either side of the B^)-T eq relation from the wind zone model 
(solid curve). This shows that the approach of using a struc- 
tural model with no ohmic heating (we refer to this as a "no 
feedback" model in ^3]) to estimate the critical magnetic field 
is a good approximation of our detailed time-evolution mod- 
els including feedback. 

Comparing the red points in Figure fl4]with the solid curve 
gives a sense of how far short the ohmicheating model falls in 
explaining the most inflated planets. For example, HAT-P-32 
is about a factor of 3-4- above the curve, so that the heating 
rate (oc B 2 ) needs to be increased by about an order of magni- 
tude to explain the observed radius. It is interesting that most 
of the unexplained objects lie within a factor of 3 in terms of 
B^o of the wind zone model. Figure 14 helps to show what 
changes to the wind zone model woula explain more of the 
observed objects. We have assumed the relation 7i so = 1.33r eq 
(from eq. ]3l)); a larger factor between 7i so and T eq would 
move the solia curve to the left, allowing ohmic heating to 
explain the radii of low T eq planets such as WASP-06. The 
dashed and dotted curves show the effect of changing B r . In- 
creasing B r from 10 to 100 G does increase B<m at low temper- 
atures, but reduces B^q at high temperatures where magnetic 
drag is enhanced. A larger depth AlnP would help to reduce 
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FIG. 13. — The ratio of observed planet radius R \, s and predicted radius 
^pred (3 Gyr) for observed hot jupiters as a function of T C q (upper panel) or 
M (lower panel). In the upper panel, the size of the circle scales with planet 
mass; in the lower panel, the size of the circle scales with r cq . 

the number of discrepant objects since both Z? a d v and fidrag in- 
crease with AlnP (eqs. (29) and (30)). 



6. SUMMARY AND DISCUSSION 

In this paper, we present models of ohmic heating in hot 
jupiters in which we attempt to decouple the interior and wind 
zone by replacing the wind zone by a boundary temperature 
7i so and magnetic field B^o, both evaluated at a pressure p = 
10 bars. This approach allows us to survey the outcomes of 
ohmic heating, parametrized by T KO and B^o for planets with 
different mass M. This is similar in spirit to models of gas 
giant cooling, which often set an outer boundary pressure of 
lObars and separately integrate a rio-r e ff relation to use as an 
outer boundary condition. 

The main conclusions of the paper are: 

1 . Figure [4] is a key result, showing as a function of en- 
tropy how the ohmic power compares to the planet luminosity. 
Only planets with entropy below a critical value have enough 
ohmic heating to slow their contraction rate. Of particular 
note are the different mass dependences: at fixed B^q and 7i so , 
the cooling luminosity L oc M whereas the ohmic power de- 
creases with mass (we find P h m oc R 2A /M). 

2. Ohmic heating has two effects on the thermal state of 
the planet. As well as providing direct heat input into the 
adiaba tic convective interior (as found by previous works, see 
|Batygin & Stevenson| ( |20T0) ; |Perna et al.| ( |2010b| l; |Wu & Lith-| 



wick (2012 1), the feedback of ohmic heating in the region 
between the wind zone and the convective boundary moves 
the convective zone boundary deeper (Fig. [6j, leading to a re- 
duced cooling luminosity and reduced internal ohmic heating. 
Because the electrical conductivity changes dramatically with 
pressure through the planet, the total ohmic power inside the 
convection zone is very sensitive to its radial extent. To com- 
puting the planet age and radius at the late stage when ohmic 
heating halted the cooling, it is crucial to accurately locate the 
convective-radiative boundary. 

3. A larger is required for ohmic heating to be im- 
portant in more massive planets or planets with larger T eq . 
This can be seen in Figures [7] and [8] which show the age of 
a cooling gas giant when ohmic heating becomes important, 
and the magnetic field strength required for ohmic heating to 
be important at different values of 7j so . For example, at a 
temperature 7j so = 1750 K, Figure [8] shows that B^q ss 30 G 
will halt the contraction of a 0.3 Mjplanet in 3 Gyr, whereas 
Bc/,0 sa 150 G is required for a IMj planet at that temperature. 
At higher 7j so = 2250 K, the required values are B^q s» 100 G 
for a 0.3M y planet or B^ ps 700 G for a IMj planet. 

4. With a specific model for the wind zone (S j5.1[ ), we can 
compare to observed systems as a function of their observed 
equilibrium temperatures r eq . The wind zone model specifies 
the induced field B$o (or equivalently, the radial current that 
penetrates into the interior; see eq. (9)) as a function of T eq , 
and the relation between r eq and the temperature at 10 bars. 
Using the scaling analysis proposed by Menou|( f2012[ ) for the 
dynamics of the wind zone, together with the atmospheric 
temperature profile from Guillot (2010 ), we find that it is diffi- 
cult for ohmic heating to explain the large radii of hot jupiters 
with large masses and large T eq (see Fig. (13)). 

5. A more general approach is to calculate, for each ob- 
served planet, the B^o that is required if ohmic heating is pro- 
viding a significant fraction of the luminosity (and therefore 
able to significantly change the contraction rate of the planet). 
This is shown in Figure 14 and shows how much the hea ting 
rate needs to be increasedover the wind zone model in §5.1| 
to explain particular objects. A modest increase in the wind 
zone thickness over that assumed here, or larger ratio of the 
temperature at depth 7j so compared to r eq , would i mpro ve the 
agreement with observed radii (see discussion in §5. 2} . Even 
so, several objects require a much more dramatic increase in 
heating rate (see Fig. 
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The difficulty in explainin g many of the observed rad ii that 
we have found differs from Batyg irTet al.| ( |20TT] l and |Wu &] 
Li th wick (2012) who found that they could account for almost 
all of the observed hot jupiter radii with ohmic heating. The 
key difference is that we do not assume here that the heat- 
ing efficiency (the fraction of the irradiation going into ohmic 
power, typically taken to be e ~ 1%) to be fixed, but instead 
use the wind zone model to set the induced magnetic field in 
the wind zone and therefore the magnitude of the heating. 

It is important to emphasize that our conclusions about the 
efficacy of ohmic heating depend on the particular prescrip- 
tion for the magnetic field in the wind zone that we have used. 
In fact, many complexities underlie the path from the irradi- 
ation to the properties of the induced magnetic field. More 
realistic 3D wind zone models may give a different picture 
than the simple ID force balance scalings we have used here. 
For example, in this paper we have assu med the wind zo ne ex- 
tends top= lObars. Figure 3 of Wu & Lithwick (2012) nicely 
illustrates the importance of the depth of the wind zone, show- 
ing that a shallower wind zone requires a significantly larger 
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FIG. 14. — For each observed planet, we show the B^q required for ohmic heating in the convection zone to be 30% of the luminosity as estimated from 
no-feedback planet models, and compare with the results of our time-dependent calculations with feedback included. The curves show the B^o-Tiq pred icted 
by the wind zone model for three different values of B r . Black p oint s show planets whose radii can be explained by our model (i? prcc | > R b s in Fig.|13|, red 
points show planets that cannot be explained (R ple d < R c bs m Fig. 1131. For clarity, we use the following abbreviations for planet names: W-WASP; Fi=HAT-P; 
K-Kepler; OG-OGLE-TR; C-CoRoT-P. 1—1 



overall efficiency to achieve the same interior heating. One 
situation in which this will break down is for young plan- 
ets with high entropies when the radiative/convective zone 
boundary is at lower pressure. More work is needed on what 
happens when the interior convection zone extends into the 
wind zone region. 

Our results emphasize the key inputs that are necessary 
from atmospheric models: the thermal structure and dynam- 
ics of the wind zone including a large scale magnetic field, 
the values of induced magnetic field, or equivalently the mag- 
netic Reynolds number Rm, that can be attained there, and the 
depth of the wind zone. More studies of the local conductiv- 
ity profile and magnetic field properties in the high magnetic 
Reynolds number regime are needed. In particular, it is not 
clear whether the large values of induced field B^q > 1000 G 
needed to explain the observed radii (Fig. 17) can be achieved 
in the wind zone. Furthermore, whether the implied large in- 
ternal currents affect the planetary dynamo is also an open 
question. 

Our results do compare favorably with previous calcula- 
tions if we use equation ( 1 1 1 to set a value of B^ appropriate 
for the wind zone conditions assumed in those papers. For ex- 
ample, we are able to compute t he 3% heating profile at pres 
sures p > 10 bars in Figure 4 of Batygi n~& Stevenson| ( |2010" 
by setting B^p = 300 G; we reproduce the heating profile of 



Tres-4b from |Wu & LithwTck| | |2012) with B 0O « 1500 G 
However, a complication in comparing different models is 
that the heating dissipated in the wind zone is coupled with 
the heating dissipated in the planet interior. Wu & Lithwick 



(2012) in particular discuss the expected ratio of heating de- 
posited in different layers. But this ratio is generally model 
dependent and varies through the planet lifetime. A direct re- 
sult of this kind of coupling is that models with same heating 
efficiency but different wind zone model are not physically 
comparable. For example , for a given set of planet properties, 
the radi us predicted b y Batygin & Stevenson] ( |2010| l is larger 
than in |Wu & Lithwick| ( |2012| l for the same choice of heat- 
ing efficiency e, because the heating ratio b etween the wind 
zone a nd the interior is much smaller in Batygin & Stevenson 
(2010), creating a muc h stronger internal heat source. Simi- 
larly, although |Menou| ( |201 2) estimated the total ohmic heat- 
ing efficiency to be > 1 % over a certain range of equilibrium 
temperatures (with the weather layer between 60 mbars and 
150 mbars), the internal heating has a much lower efficiency, 
consistent with our findings in ^5] 

Another uncertainty is in the microphysics aspects of the 
electrical conductivity. For example, as we noted in §3.1, 
Batygin & Stevenson (2010) make a different choice for the 
electron-neutral cross-section and thermal averaging that re- 
sults in a factor of 9 difference in electr ical conductivity than 
we adopt here. The estimates in 5 2.2 show that the amount 
of ohmic power is sensitive to changes in the electrical con- 
ductivity (or the ionization fraction) in two ways. At low den- 
sities in the wind zone, the conductivity determines the size 
of the current (eq. flO)); in the interior, the ohmic power is 
oc l/rj (eq. 1 13 1). For a fixed efficiency e, a different nor- 
malization for a does not change the evolution of the planet, 
since the normalization of the heating profile is determined 
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FIG. 15. — The opacity profile in a planet with parameters as in Table 
1 tablenote a (no ohmic heatin g, radiative model) , sh owing the opacity as 
calculated by combining the Freedman et al. 1 2008 1 and Potekhin & Chabrier 
H |2010> tables (solid curve) or trom the MESA code ^Faxton et al.|2UI I) (red 
dashed curve). 
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FIG. 16. — The contribution to the electron fraction Y e from different alkali 
metals as a function of temperature and pressure. Solid curves are for potas- 
sium, dotted for sodium, and dashed for aluminum. In each case, we show 
(top to bottom) pressures of 1, 100 and 1000 bars. 



by the choice of e, and a(r) determines only its shape. The 
normalization of a is important, however, when going beyond 
the constant efficiency assumption, making it crucial to under- 
stand the processes that set the ionization level in hot jupiter 
atmospheres. 



This work began as a project at the 201 1 International Sum- 
mer Institute in Modeling in Astrophysics (ISIMA), held at 
the Kavli Institute of Astronomy and Astrophysics, Beijing, 
China. We thank ISIMA for support and KIAA for hospital- 
ity during the program. We would like to thank E. Chiang, 
D.N.C. Lin, A. P. Showman, Y. Wu and Y. Lithwick for use- 
ful discussions during 2011 ISIMA. We are also grateful for 
the helpful suggestions during private communication from 
T. Guillot and R. Laine, and to G.-D. Marleau for discussions 
about gas giant models and a thorough reading of the paper. 



10" 



10"' io"' io" 10 io' io J io' io° 10 10' 

Pressure p (bar) 



10 




10 10' io° 10 

Pressure p (bar) 

FIG. 17. — Top panel: the electrical conductivity profile of a planet with 
parameters as in Table 1 tablenote a (no ohmic heating, radiative model), 
showing the contributions from alkali metals (dashed green curve) and hy- 
drogen (dashed blue curve). Bottom panel: the electron fraction Y e as a func- 
tion of pressure, with the contribution from alkali metal ionization shown as 
a dashed curve. 
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APPENDIX 

MICROPHYSICS OF THE PLANET INTERIOR 

We discuss the microphysics input in our gas giant mod- 
els her e. We adopt the equation of state from Saumon et al. 
( |1995[ ) with helium fraction Y = 0.25. In order to maximize 
the planet radius, we do not include a solid core or elements 
heavier than helium. 



The radiative opacity is taken from Freedman et al. ( 2008 1 
and in the core we include therm al conduction by electrons 
from |Potekhin & Cha brier (2010). The transition from radia- 
tive to conducive opacity occurs at a pressure which is greater 
than the maximu m pressure of 300 bars covered by the Freed- 
man et al. ( 2008[ ) tables. In the intermediate regime, we as- 
sume the scaling k oc p 5 . The opacity profile for our stan- 
dard model is shown in Figure [T3J over-plotted wi th opacity 
taken f rom MESA using the same planet structure (Paxton et 
|al.|201ip . 

The electrical conductivity has contributions from alkali 
metal ionization in the outer layers, and hydrogen in the inte- 
rior. In the upper atmosphere of hot jupiters, the conductivity 
is set by the ionization of alkali metals. For potassium, which 
has the lowest ionization potential^] The potassium only Saha 

5 The first ionization potentials of K, Na, Al, Mg and Fe are 4.34, 5.14, 
5.99, 7.65 and 7.90 eV respectively (David 2003]. 
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equation (Balb us~& Hawley||2000 Perna et al. 2010a) gives 
the ionization fraction xt = n e /n as 



x k = 



m e k B T \ 
2-kTl 1 J 



3/2 



-4.35eV/k B T 



(Al) 



= i.03xio- 3 r 3 5/ V 2519 ^ 



Jk_ 

10- 7 



1/2 



(—) 

llbar/ 



-1/2 



where f K is the number fraction of potassium. Although 
potassium dominates, we also include the contribution of Na, 
Mg, and Fe in the ionization balance to sum up the total ion- 
ization fraction. The ionization fraction of each alkali metal 
is computed separately by assuming a balance independent 
on the presents of other elements . We do not include the con- 
tribution of Al in t he calculation, which is likely condensed 
out ( Lodders|1999| l. But our results are not very sensitive to 
elements beyond potassium. This is illustrated in Figure 16 
which shows the contributions to the ionization level from 
Na and Al at different pressures. Once the ionization fraction 
is determined, the conductivity is a = n e e 2 /m e v where the col- 
lision freq uency of electron-n eutral collisions is v w = n n (av) e 
given by |Draine et al.| (fT983) as 



In the deeper part of the planet, the hydrogen is ion- 
ized by high pressure and the conductivity is dominated by 
electron-proton collisions. In the fully-degenerate limit, v ep i = 

4e 4 m e A/3irft 2 = 1.8 x 10 16 s _I . In the non-degenerate limit, 
t'epnd = 6.4 x 10 23 s' 1 pY e T~ 3 / 2 , in which Y e is the electron frac- 
tion. We interpolate between the two limits to give an esti- 
mation of the total contribution: v~ 2 = v~ 2 & + v~ 2 & . We also 
include the cond uctivity at intermediate densities as given by 
|Liu et al.] p006 ). Before the hydrogen molecule is fully ion- 
ized, the band-gap of hydrogen will diminish with increasing 
pressure, to a lev el where there is a significant contribution to 
the conductivity. Liu et al.| ( [200 6) give this as 



cr s = <7 exp 



-E g (p) 
k B T 



(A4) 



(av) e = 10" 



. 1S ( \2%k B T 
\ 9irm e 



1/2 



3 -1 
cm s . 



(A2) 



The conductivity is then 



1500 K 



-1/2 



(A3) 
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where between 0.2 Mbars and 1.8 Mbars, E g = 20.3-64.7p, 
where E g is in eV, and p is in molcm -3 , and 00 = 3.4 x 
10 20 s" 1 exp(-44p). The overall conductivity is constructed 
as a = a s +n e e 2 /m e v = a s + 1.52 x l(P 2 pY e /i/. The collisional 
frequency v is the sum of electron-neutral and electron-proton 
collisions. A typical conductivity profile and the contribution 
of different components are shown in Figure [TT] 
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